% ! problem9.m --> based on <<MATLAB codes for Finite Element Analysis>>

clear ;
clc
format default
p=struct();
Case=3;

% define model parameters
% unit:SI(mm)
p.EI = 1;
p.L = 1;
% * distribute force
p.P = -1;
% define elems num of bornoulli beam
p.elem_num=80;
% define problem's dimension : 1D/2D/3D;
%        = each node's dof  
p.problem_dimension=2;

% ? -------------------------------

% node's coordinates (x)
p.nodes=linspace(0,p.L,p.elem_num+1)';
p.elems=zeros(p.elem_num,2);
for i=1:p.elem_num
    p.elems(i,1)=i;
    p.elems(i,2)=i+1;
end
p.node_num=size(p.nodes,1);
p.node_Coord_x=p.nodes(:,1);

% global degree of freedom number
p.global_dof_num=2*p.node_num;

% define global Nodal displacement  colum vector
p.displacements=zeros(p.global_dof_num,1);
p.node_forces=zeros(p.global_dof_num,1);

% define boundary conditions
if Case==1
    % Case 01: clamped at x=0
    % fixed dof 
    p.fix_dof=[1 2]';
end
if Case==2
    % Case 02: clamped-clamped
    % fixed dof 
    p.fix_dof=[1 p.elem_num*2+1 2 p.elem_num*2+2]';
end
if Case==3
    % Case 03: simply supported-simply supported 
    p.fix_dof=[1 2*p.elem_num+1]';
end
% p.prescribed_dof=sort([p.fix_dof;p.fix_dof]); % 排序

% initial global stiffness matrix
p.global_stiffness_matrix=zeros(p.global_dof_num);

% compute all ElemStiffnessMatrix and assembly stiff matrix
elem_global_dof_num=p.problem_dimension*size(p.elems,2);
p.elemStiffs=zeros(elem_global_dof_num,elem_global_dof_num,p.elem_num);

for i=1:p.elem_num
    connectivity=p.elems(i,:);
    % elem's all dof : a 3d truss node has 2 node,each node have 2 dof,
    % i-th node --> w_(2*i-1) / dw_(2*i-2) /  dof(displacement of node) 
    elem_dof=[connectivity(1)*2-1 connectivity(1)*2  ...
        connectivity(2)*2-1 connectivity(2)*2 ];

    node1_x=p.node_Coord_x(connectivity(1));
    node2_x=p.node_Coord_x(connectivity(2)); 
    
    elem_length=node2_x-node1_x;
    a=0.5*elem_length;
    k_e=(p.EI/(a^3))*0.5*[3 3*a -3 3*a;
                        3*a 4*a*a -3*a 2*a*a;
                        -3 -3*a 3 -3*a;
                        3*a 2*a*a -3*a 4*a*a];
    p.elemStiffs(:,:,i)=k_e;
    p.node_forces(elem_dof)=p.node_forces(elem_dof)+((p.P*a/3)*[3 a 3 -a])';

    % assemble stiffness matrix
    p.global_stiffness_matrix(elem_dof,elem_dof)=...
        p.global_stiffness_matrix(elem_dof,elem_dof)+k_e;
end

% apply boundary condition 
p.displacements(p.fix_dof)=0;
% p.node_forces(p.load_dof)=p.P;

% solution KU=F
p = solutionStruct(p);

plot(p.node_Coord_x,p.displacements(1:2:2*p.node_num),'.')

minDispalce=min(p.displacements(1:2:2*p.node_num));
text(0.35,0.5*minDispalce,"挠度最小值:"+minDispalce)

if Case==1
    title("Case 01: clamped at x=0")
end
if Case==2
    title("Case 02: clamped-clamped")
end
if Case==3
    title("Case 03: simply supported-simply supported") 
end


% function new_p=solutionStruct_bornoulli_beam(p)
%     % free dof in structure
%     free_dof=setdiff((1:p.global_dof_num)',p.fix_dof);
    
%     % solute the equation :KU=F --> U=K\F
%     U=p.global_stiffness_matrix(free_dof,free_dof)\p.node_forces(free_dof);
    
%     p.displacements(free_dof)=U;
%     % displacements=zeros(p.global_dof_num,1);
    
%     p.node_forces_fu=p.global_stiffness_matrix*p.displacements;
    
%     new_p=p;
% end


















